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The effects of particle loading on 
turbulence structure and modelling 

By K. D. SQUIRES AND J. K. EATON 

1. Introduction 

Conventional wisdom implies that the presence of particles provides an ad- 
ditional dissipation or attenuation of turbulence. However, it is not clear how 
this extra dissipation may be incorporated into turbulence models or how it de- 
pends on the parameters of the problem including the particle time constant, 
the mass loading, and the various dimensionless parameters describing the tur- 
bulence. Simply focusing on the extra dissipation completely neglects the effect 
of turbulence structure on the instantaneous particle concentration field and 
the possibility of interactions between the particle cloud and instability mecha- 
nisms generating turbulence. For example, particles selectively concentrated in 
particular structures may cause rapid attenuation of that structure or trigger a 
new instability mechanism. Therefore, it is important to determine the effect of 
turbulence structure on the behavior of the particle cloud. 

The objective of the present research was to extend the DNS approach to 
particle-laden turbulent flows using a simple model of particle/flow interaction. 
The program addressed the simplest type of flow, homogeneous, isotropic tur- 
bulence, and examined interactions between the particles and gas phase tur- 
bulence. The specific range of problems examined include those in which the 
particle is much smaller than the smallest length scales of the turbulence yet 
heavy enough to slip relative to the flow. The particle mass loading is large 
enough to have a significant impact on the turbulence, while the volume load- 
ing was small enough such that particle-particle interactions could be neglected. 
Therefore, these simulations are relevant to practical problems involving small, 
dense particles conveyed by turbulent gas flows at moderate loadings. 

This report presents a sample of the results illustrating modifications of the 
particle concentration field caused by the turbulence structure and also illus- 
trating attenuation of turbulence by the particle cloud. 

2. Overview of the simulations 

The numerical method used in the present research is based on the pseudo- 
spectral technique developed by Rogallo (1981) for solving the incompressible 
Navier-Stokes equations. A modification of Rogallo’s code was used to simulta- 
neously track a large number of particles using a vectorized, second-order inter- 
polation scheme. This code is time advanced using second-order Runge-Kutta, 
and at each substep the fluid velocity at the particle position was calculated 
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using trilinear interpolation. Numerical experiments showed that more accurate 
interpolation schemes do not significantly improve the results. The new position 
was obtained by then integrating the particle equation of motion 


dvj(t) _ Uj(Xj(t),t) - Vi (t) 
dt t p 
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where X 3 (t) and v,(t) are the position and velocity of the particle, r p is the 
particle time constant, and Ui(X 3 (t),t) is the velocity of the fluid. Stokes’ law has 
been used to calculate the drag on a particle, and it is, therefore, assumed that 
the particle is much smaller than any lengthscale of the flow (the Kolmogorov 
scale, rj). It is also assumed that the particles are sufficiently dense such that 
other forces in the full equation of motion, e.g. buoyancy and added mass, are 
negligible compared to the Stokes drag. Finally, the particles are assumed to 
occupy a negligible volume fraction, and particle-particle interactions are also 
assumed negligible. 

To calculate the effect of the particles on the turbulence, the momentum 
equation of the fluid was modified by a source term on the right-hand side 
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where c(x,-, t) is the particle mass per unit of volume, or particle density (repeated 
indices imply summation). 

The problem studied is homogeneous, isotropic turbulence in a box. “Natural” 
isotropic turbulence decays, and the lack of stationarity complicates analysis of 
the results. Therefore, the low wave number modes (large scales) are artificially 
forced to maintain stationarity using the scheme developed by Hunt, et al (1987). 
Results were obtained using 373000 particles for simulations using 32 s points and 
l.xlO 6 particles for simulations using 64* points. The relatively coarse grid 32 s 
calculations permitted a parameter study in particle time constant and mass 
loading ratio. For the 32* simulations the values of the particle time constant 
divided by the eddy turnover time used were 0.14, 0.75, and 1.50. The mass- 
loadings used in the 32* simulations were 0.1, 0.5, and 1.0. All three values of 
the mass loading were used for each of the time constants for these simulations. 
For the 64* simulations the values of the dimensionless time constant were 0.075, 
0.15, and 0.52 and the mass loadings were 0.1, 0.5, and 1.0. 

For each simulation the two phases were uncoupled and time-advanced to reach 
a statistically stationary state and eliminate initial transients of the forcing. At 
this point in time, the momentum equation of the fluid was coupled to the 
particle momentum equation with a specified mass loading. Once the phases 
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were coupled, the simulations were time advanced just over 12 eddy turnover 
times for the 32 s simulations and slightly over 6 eddy turnover times for the 64 s 
simulations. 

Results from the simulations of no coupling (4> = 0 where <j> is the mass loading) 
between the phases provide a baseline case to compare the effect of increased 
mass loading on the turbulence statistics. These simulations are also useful 
for examining the effect of turbulence “structure” on the particle concentration 
field. For these baseline simulations the Reynolds number based on the Taylor 
microscale was approximately 37. This relatively low value for the Reynolds 
number assured good resolution of the velocity field. 

3. Results and discussion 

Maxey (1987) has shown that particles will tend to accumulate in regions of 
high strain rate or low vorticity. This effect can be measured by using the second 
invariant of the deformation tensor, dui/dxj. The second invariant of dui/dxj 
is given by 

I Id = — ~(5 a — — Wj-Wy) (3-1) 

where s 1 is the magnitude of the strain rate, , u»y is the jth component of 
the vorticity vector and wyu/y is by definition the enstrophy. 

Thus, from equation (3.1) highly vortical regions correspond to I Id being large 
and positive (where the number density is small), while regions of high strain 
correspond to I Id being large and negative (where the number density is large). 

Contours of number density for each of the particles used in the simulations 
showed distinct regions where particles had collected. It was found that the 
peak number density is as much as 40-50 times the mean value for the lightest 
particles used in the simulations. By comparing these contours with the contours 
of the second invariant of dui/dxj , it was found that the particles showed a 
tendency to accumulate in regions of low vorticity or high strain rate. As the 
particles were made more sluggish (i.e., larger r p ), there is less of a tendency 
for them to accumulate in these regions. For the heaviest particle used in the 
32* simulations (r p /T e = 1.50), the particle concentration field was found to be 
essentially random. 

The preceding results were quantified by examining the conditional expecta- 
tion of the number density given the value of the enstrophy. The conditional 
expectation showed that the effect of increasing particle inertia reduced the like- 
lihood of finding the particles in regions of low vorticity (figure 1). For the 
64 s simulations, the lightest particles (r p /T e = 0.075) actually showed less of 
a tendency to accumulate in regions of low vorticity than did particles which 
were twice as heavy ( r p /T e — 0.150). This can be explained by the fact that 
very light particles can follow nearly all of the turbulent motions and would, 
therefore, show no preference to be in regions of low vorticity. 
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<p — 0.0 



FIGURE 1. Conditional expectation of number density given enstrophy for 
the 64* simulations, <j> = 0. 


For increasing values of the mass loading, the conditional expectation of num- 
ber density given the enstrophy for the light particles showed that they had less 
of a tendency to be in regions of low vorticity. However, this was not the case 
for the heavier particles. 

Correlations between the number density and enstrophy showed that as the 
loading was increased the lightest particles became less correlated with enstro- 
phy. For the larger time constants, however, the number density and enstrophy 
become more correlated for increasing values of the time constant and increased 
loading. Correlations between number density and II<x showed that as the mass 
loading was increased the correlation between number density and lid increased. 
This increase in the correlation was more significant for larger values of the time 
constant. 

Once the two phases were coupled, the turbulence evolved to a new stationary 
state within about two eddy turnover times. For all the time constants used in 
the simulation, it was found that the time required to come to a new equilibrium 
was longer for increasing loadings, and the maximum development time for any 
case was just over two eddy turnover times. Time averaged statistics such as 
turbulence energy showed that increased mass loading decreased the turbulence 
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FIGURE 2. Turbulence energy versus loading from the 32 s simulations, 
r p /T e = 0.75. 

energy (figure 2). These results are consistent with the fact that the drag of the 
particles on the turbulence acts as an additional source of dissipation. This can 
be shown be deriving the transport equations for (beginning with equation 
(2.3)). 

Frequency spectra measured along the particle path were nearly identical for 
all simulations. The differences in the turbulence spectrum for the three cases 
can be attributed to the fact that the lighter particles are found more often in 
regions of low vorticity and high strain rate than are the heavier particles. This 
causes a sampling bias resulting in the differences between the spectra along the 
particle path. Measured values of the particle velocity spectra were found to be 
in excellent agreement with a theoretical prediction of Csanady (1963). It was 
somewhat surprising to find such close agreement with Csanady’s theory in view 
of the fact that the particles are concentrated within particular regions of the 
turbulence. Examination of the power spectra of the turbulence for increasing 
values of the mass loading showed that energy was removed nearly equally from 
all frequencies. This may be due in part to the low Reynolds number of the 
simulations. Selective energy loss within certain frequency bands may occur at 
higher Reynolds numbers. 
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Correlations between enstrophy and pressure were decreased more by the light 
particles than by the heavy particles for increasing values of the mass loading. 
Since regions of low pressure are associated with regions of high vorticity, these 
results indicate that the lighter particles cause more distortion of these eddies 
than do the heavier particles. The fact that the lighter particles show a more 
pronounced tendency to accumulate in regions of lower vorticity (for tfi = 0) than 
do the heavier particles may provide some explanation for this. 

Correlations between I Id and pressure showed that for the time constants and 
mass loadings used in the simulations the correlation between these two quanti- 
ties remained reasonably constant, decreasing slightly for the lightest particles. 
Correlating lid with pressure correlates regions of low pressure with vortices 
(large positive lid) and regions of high pressure with straining regions (large 
negative I Id)- Therefore, these results indicate that the correlation between the 
straining regions with regions of high pressure must be increasing to compen- 
sate for the loss of correlation between the vortical regions and regions of low 
pressure for the lighter particles. 

In summary, it was found that for the case of zero loading there are significant 
effects of the turbulence “structure” on the resulting concentration field. These 
results were quantified by measuring conditional expectations and correlations. 
It is shown that the lighter particles show a strong tendency to be in regions of 
low vorticity and high strain rate. 

For increasing values of the mass loading the power spectra of the turbulent 
fluctuations do not show any frequencies to be preferentially modified by the 
particles. This may be partially due to the fact that there is not a large range 
of scales in the simulations. 

The lightest particles used in these simulations were found to modify the 
turbulence field differently than did the heavier particles. Since the light particles 
show a more pronounced tendency to accumulate in regions of low vorticity for 
<j> = 0, these regions are modified more by the light particles as the loading is 
increased than the heavy particles. Evidently, this selective modification by the 
light particles causes more of a distortion of the turbulent eddies than the more 
uniform modification by the heavier particles. 
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